Interaction-driven equilibrium and statistical laws in small systems. The cerium atom. 
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I. INTRODUCTION 
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ON ' It is shown that statistical mechanics is applicable to isolated quantum systems with finite 

ON 1 numbers of particles, such as complex atoms, atomic clusters, or quantum dots in solids, where 

the residual two-body interaction is sufficiently strong. This interaction mixes the unperturbed 
^ ' shell-model (Hartree-Fock) basis states and produces chaotic many-body eigenstates. As a result, 

an interaction-induced statistical equilibrium emerges in the system. This equilibrium is due to the 
off-diagonal matrix elements of the Hamiltonian. We show that it can be described by means of 
temperature introduced through the canonical-type distribution. However, the interaction between 
the particles can lead to prominent deviations of the equilibrium distribution of the occupation 
numbers from the Fermi-Dirac shape. Besides that, the off-diagonal part of the Hamiltonian gives 
rise to the increase of the effective temperature of the system (statistical effect of the interaction). 
For example, this takes place in the cerium atom which has four valence electrons and which is used 
in our work to compare the theory with realistic numerical calculations. 
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Consider a many-body quantum system of interacting particles. If the number of particles is large, statistical laws 
can be applied to describe the properties of the system. They can also be applied to few-particle systems (or even 
single particles) interacting with a heat bath. In both cases the equilibrium is achieved at arbitrarily weak interaction 
between the particles, or with the heat bath. If the interaction between the particles is strong enough, a different kind 
of statistical equilibrium is possible in isolated few-particle quantum systems. It is achieved when the excited states 
of the system become chaotic compound states. The systems examined so far are the rare-earth atom of Ce Jl) with 
just four active valence electrons, 12 nucleons in the s — d shell and n — 4-7 particles interacting by means of a 
two-body random interaction . In spite of the obvious differences these systems have much in common, as far as 
properties of their eigenstates and various statistics are concerned. It has been shown in that in the regime of 
' compound excited states one can introduce such thermodynamic parameters as temperature and entropy, and observe 
£NJ , other typically statistical features, e.g., the Fermi-Dirac distribution of the occupation numbers. 
t— I ■ In the present work we concentrate on the statistics of the occupation numbers in a realistic Fermi system: the atom 
of Ce. We show that when the interaction between the particles is strong and the two-body matrix element fluctuates 
strongly as function of the single-particle states involved (hence, there is no good mean-field approximation), large 
deviations from the Fermi-Dirac behaviour are observed. However, a statistical description of the system including 
the introduction of a temperature is still possible if the interaction between the particles is properly accounted for. 

The notion of compound states is important for our work, so we would like to explain it in greater detail. Suppose 
that for a given range of excitation energies the interaction between the particles gives rise to a certain mean field. 
'-^J ■ This mean field can then be used to generate a set of single-particle orbitals. The multiparticle basis states of the 
system can be constructed from these orbitals by simply specifying their occupation numbers. The spectrum of such 
states in a system with several active particles is very dense since there are many ways of distributing them among 
the orbitals. For the interacting particles these multiparticle states are not the eigenstates of the system. Instead, 
they are mixed together by the residual two-body interaction. If this interaction is strong, the number of basis states 
"involved" in almost every eigenstate of the system becomes very large (about 100 in atoms and up to 10 6 in nuclei), 
and the eigenstates become almost random (chaotic) superpositions of the basis states, devoid of any good quantum 
numbers, save the exact ones - energy, parity, total angular momentum, etc. Following the nuclear physics terminology 
we call such eigenstates compound states. Their statistical properties, e.g., distribution over the unperturbed basis 
states, are very similar in different systems studied: atoms, nuclei, or a two-body random interaction model. Most 
importantly, they provide a good starting point for developing a statistical theory for isolated few-particle systems 
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II. STATISTICS OF THE OCCUPATION NUMBERS 



Statistical behaviour is usually established in the limit of a large number of particles n. Moreover, simple quantita- 
tive results can be obtained if correlations between the particles are somehow weak. This means that the interaction 
between the particles can be neglected, or - more realistically - an appropriate mean field theory is developed. The 
latter results in the picture of free quasiparticles moving in the effective self-consistent field created by the constituents. 

In the limit of large n the temperature T is a well-defined physical quantity and all equilibrium characteristics can 
be found by applying the canonical (Gibbs) probability distribution Wi oc exp(— J?W/T), where E^ 1 ' is the energy 
of the ith eigenstate of the system. For example, for a gas of noninteracting fermions this results in the famous 
Fermi-Dirac distribution (FDD) of the occupation numbers: 

(i) 



exp[( £ct - M )/T] + 1 



where e a is the energy of a single-particle state a, and /i = /i(n, T) is the chemical potential. It is determined, at a 
given temperature, from the normalization condition ^2 a n a = n. Formula ([!]), when it is valid, in fact provides one 
with a relation between the temperature and the energy E of the system: e a n a — E, which can also be viewed 
as a possible definition of the temperature. Note that Eq. (|l|) can hold for the interacting fermions as well, provided 
we consider the distribution of quasiparticles, and replace e a with the quasiparticle energy e a , which in turn depends 
on the distribution of excited quasiparticles. This is an important point of Landau's theory of Fermi liquids (see e.g., 
!)■ _ 

Strictly speaking, the FDD is derived for the grand canonical ensemble, where /x is fixed, and n is the mean number 
of particles H. For large n the difference is negligible unless fluctuations of the number of particles are concerned. 
In appendbcjA] we consider the approximations one has to make to arrive at the FDD (Q) for a finite system of n 
interacting particles obeying the canonical distribution. 

In reality there are many complex systems, such as compound nuclei, rare-earth atoms, molecules, atomic clusters, 
or quantum dots, which do not satisfy the conditions for Eq. (|l|) to hold. However, their complexity suggests that 
some statistical methods can be developed, and in nuclear physics such statistical temperature-based description has 
been known for quite a while. Intuitively such description is very natural, and a more rigorous justification does not 
seem to have been needed. 

The number of active particles in these systems can be relatively small (j$ 10), whereas the interaction between 
them (even the residual interaction in the mean-field basis) is large, i.e., greater than the energy intervals between 
unperturbed many-particle basis states. This interaction makes up for the absence of a heat bath, and promotes 
the onset of "randomization" and quantum chaos. This chaos is purely dynamical, in the sense that the Hamilto- 
nian matrix of the system does not contain any random parameters, yet, the behaviour can be complicated enough 
("chaotic"), and a number of properties, e.g., the energy level statistics, are consistent with the predictions of random 
matrix theories ■ This gives one a possibility to talk about some kind of equilibrium in the system, and pursue 
the development of a statistical theory for few-body Fermi systems jf|. 

In what follows we will look at the results obtained numerically in a realistic model of the Ce atom which contains 
only four active particles (valence electrons). We will see that the energy dependence of the occupation numbers 
differs prominently from what one expects from the FDD ([!]), and show that this behaviour results from the strong 
fluctuations of the two-body Coulomb interaction for different orbitals. We then show that this interaction can be 
taken into account within the statistical approach to calculation of the occupation numbers and other mean values, 
leading to a good agreement between the results of our statistical theory and the numerical calculations. Such 
agreement confirms the existence of equilibrium similar to the thermal one in the system of a few strongly interacting 
particles. 



III. THE CE ATOM 



The cerium atom has one of the most complicated spectra in the periodic table. The density of energy levels with 
a given total angular momentum and parity J* reaches hundreds of levels per eV at excitation energies of just a few 
eV, well below the ionization threshold of I = 5.539 eV |0,|J- The Ce atom has four valence electrons, and a well 
defined Af6s 2 5d ground state. However, with the increase of the excitation energy and involvement of yet another 
low- lying electron orbital 6p, the atomic eigenstates become compound states (in the sense of Sec. |j|), and it becomes 
absolutely impossible to choose any reasonable coupling scheme or provide any classification for them. At the same 
time, the orbital occupation numbers move away from integer values, and even the idea of a dominant configuration 
for a particular energy level becomes meaningless 0] . 
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In the present work we continue to study the cerium atom numerically, with emphasis on the energy dependence of 
the occupation numbers. The electronic structure of the Ce atom consists of a Xe-like Is 2 . . . 5p 6 core and four valence 
electrons. A large difference in the energy scales of the core and valence electrons allows us to neglect excitations from 
the core and consider the wave function of the core as a "vacuum" state |0). Accordingly, the four active electrons 
added to this vacuum form the spectrum of Ce at excitation energies E < / M . 

The calculations are performed using the Hartree-Fock-Dirac (HFD) and configuration interaction (CI) methods 
(see [jj] for details). A self-consistent HFD calculation of the neutral atom results in the construction of the core and 
valence orbitals. It also determines the mean-field potential, which is then used to calculate the basis set of single- 
particle ortho- normalized relativistic states \a) = \nljj z ) with energies e a . This procedure defines the zeroth-order 
Hamiltonian of the system, 

m =^2e a ala a . (2) 

a 

The unperturbed multi-particle basis states \k) constructed from the single-particle states, \k) — a^a^aj, o^ 4 |0), 
are eigenstates of H^: H^\k) — E^\k), where 

*f = (3) 

a 

is the zeroth-order energy of the state |fc), and ria — (k\n a \k) — (fc|aJj,a Q |fc) are the occupation numbers equal to 
or 1, depending on whether the state a is occupied in \k), or not. To subtract additional symmetries only the basis 
states |fc) with a given projection of the total angular momentum J z and parity are considered. 

The total Hamiltonian H of the active electrons is the sum of the zeroth-order mean-field Hamiltonian and 
the 2-body residual interaction 

V = - £ . (4) 

The residual interaction V contributes to the diagonal and off-diagonal matrix elements between the multiparticle 
states |fc). Its diagonal part shifts the energy of the basis state k, 

AE k =V kk = J2 (Yafipa - V aPaP ) n^nf . (5) 

a>/3 

The off-diagonal matrix elements V k i k = (^'1^1^) are responsible for mixing of the multiparticle basis states. 

Complete diagonalization of the operator H = Hd + V in the space of the basis states |fc) yields "exact" energies 
and stationary states \i), 

H\i) = E^\i), (6) 
which can be presented as superpositions of the unperturbed basis states, 



In this work we included 14 relativistic subshells nlj in the single-particle basis (6s, 7s, 6p, 7p, 5d, 6d, 4/, and 5/), 
and performed exact diagonalization of the NxN Hamiltonian matrix in a Hubert space with N ~ 8 x 10 3 , obtained 
by truncating the complete set of the shell-model atomic configurations. Our numerical results are obtained for the 
even states of Ce with the total angular momentum projection set to J z — 0. Thereby, all possible states with J from 
to 10 are included. For the given choice of the basis the numbers of eigenstates with J = 0-10 are 343, 917, 1354, 
1493, 1433, 1153, 826, 497, 262, 107, and 34, respectively. 

The eigenvalue densities pj(E) for J = 4-8 are shown in Fig. |l|. They have been window-averaged over AE = 0.05 
a.u. to smooth out short-range fluctuations. The largest density is observed for J = 4, and with the exception of small 
regions near the ends of the spectra, all pj(E) are proportional to each other. The shapes of the eigenvalue densities 
are basically determined by the corresponding basis state densities (although the effect of level repulsion makes the 
former slightly wider). They are characterized by a very rapid increase in the low-energy part. This increase is a 
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direct consequence of the fact that the accessible energy can be distributed in an ever greater number of ways between 
the 4 electrons. Being essentially of combinatorial nature, the level density can be described by the exponent 



p{E) cx exp [ay/E - E g j (8) 

which is derived in the independent-particle model Jl(J [E g is the ground state energy). As seen from Fig. [j] the level 
density from our calculation indeed follows (H) at low energies, but then reaches its maximum and decreases [ jTl| . 
These latter features arc unphysical as they are consequences of the finite size of our basis. However, within 5 eV of 
the ground state the configuration set we use is reasonably complete. 

Relativistic atomic subshells nlj are (2j + l)-degenerate, therefore, we consider average occupation numbers 

n s =97 1 ^ ^ 9s 1 ^2ala a , (9) 

aifunctionns qGs 

where g s — 2j + 1 is the degeneracy of the subshell s. In the system of a large number of weakly interacting particles 
thermally averaged values of n s are given by the FDD (|l|) . In the quantum dynamical system, like the Ce atom, the 
occupation numbers for any eigenstate can be obtained as ni = (i\n s \i). When the number of active particles is 
small n s show strong lcvcl-to-lcvcl fluctuations, and it is more instructive to look at the spectrally averaged values 



n s (E) = {i\n,\i) = \ C k } ( fc l"sl fc > . ( 10 ) 

k 

where overline means averaging over the eigenstates i within some energy interval around E. 

A typical distribution of the occupation numbers calculated at the excitation energy of 3.75 eV above the atomic 
ground state is shown in Fig. || as a function of the single-particle energy e s of the orbitals (see Table |). The values 
of e s are determined with respect to the Xe-like Ce 4+ core. One can see that the distribution does not look at all like 
a monotonically decreasing FDD. A similar picture is observed over the whole energy interval from the ground state 
to the ionization potential. For example, the lowest even state of Ce has a configuration of 4/ 2 6s 2 , while the FDD 
would tell us that all 4 electrons must be placed in the lowest 4/-orbital, when the energy or "temperature" of the 
system is low. In reality the 4/ 4 electron configuration lies at very high energies due to a strong electron repulsion in 
this compact orbital (the radius of the 4/ orbital in Ce is at least two times smaller than that of any other valence 
orbital). 

Of course, considering the orbital energies e s has the drawback that they completely ignore the residual interaction 
between the valence electrons. When this interaction is strong one would wish to introduce some new mean-field orbital 
energies e s that would incorporate the effect of such interaction. The value of e s for the orbital s will inevitably depend 
on the distribution of the other electrons, and hence on the excitation energy of the system. In Sec. [v| we introduce 
such energies within the statistical approach. However, even when the occupation numbers are plotted against e s 
there is still a large deviation from the standard FDD distribution. 

At first sight such a strong deviation from the FDD in a strongly interacting Fermi system speaks against any 
possibility of a statistical description of the system. In what follows we show that the strong Coulomb interaction 
between the electrons can be incorporated in the canonical-ensemble description of the system, and thermally averaged 
occupation numbers can be n s (T) derived. As a result of the electron interaction, the distribution of n s (T) differs 
from the FDD. What is more important, we find good agreement between the results obtained by means of this 
statistical approach and those from the pure dynamic calculation of the Ce eigenstates, n s (T) « n s (E). We also find 
a way to relate the energy and the temperature. 



IV. CANONICAL ENSEMBLE AND THE STRENGTH FUNCTION 



In this section we show that averaging over the canonical distribution, which weighs different states according to 
their energies E^ with probabilities Wk <x exp(— Ek/T), is very similar to averaging over the exact eigenstates \i), 
when these eigenstates are compound, i.e., include large numbers N c of basis states \k) mixed together with small 

weights Cfc ~ 1/y/Nc, Eq. (Q). For a classical system the latter is equivalent to averaging over the microcanonical 
distribution that considers all points on the surface E — const of the phase space as equally probable p"2| . As is known, 
the two types of averaging yield identical results for large systems [|5| . 
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A. Averaging over the canonical distribution at a given temperature. 



Suppose first that the off-diagonal part of the residual two-body interaction V is switched off. The multiparticle 

basis states |fc) then correspond to the stationary states of the system with energies Ek = + AE^. The interaction 
with a heat bath at temperature T results in the canonical distribution of probabilities of finding the system in a 
given state k, Wk = Z^ 1 exp(— Ek/T), where Z = J^ fc exp(— Ek/T), so that J2k Wk = 1- The occupation numbers at 
temperature T are calculated as 

n s (T) = J2w k (k\n s \k) = Z" 1 ^exp(-S fc /T)(fc|ri fl |A:) . (11) 

k k 

The spectrum of Ek is similar to the eigenvalue spectra shown in Fig. [I], and is characterized by a rapid rise of its 
density [see Eq. (0)]. Thus, if we replace summation in Eq. (O) by integration over Ek, 



n s (T) = I w T (Ek)(k\h s \k)dE k , (12) 

where we introduced the probability density WT(Ek) = Z^ 1 exp(—Ek/T)p(Ek), the integrand in ( |l2| ) will peak strongly 
due to competition between the two exponents, the decreasing exp(— Ek/T) and the rising p{Ek). As a result, the 
main contribution to n s (T) is given by the vicinity of the maximum of wt (Ek ) ■ The equation for the position of the 
maximum Ek = E provides a relation between the most probable energy E and the temperature, 

1 d{\n\p(E)}} , , 

-_ + Aigp =0 . (13) 

If the temperature is not too small the maximum of WT(Ek) is almost symmetric, and the most probable energy 
becomes close to the mean energy: 



E » E(T) = J E k w T {Ek)dE k . (14) 
If we use the analytical form (||), Eq. (|l3|) yields 

T=- a ^E~~E g . (15) 



B. Strength function and averaging over the compound states. 



Let us now come back to the dynamic description of the isolated many-body quantum system and switch on the 
off-diagonal part of the residual interaction. In this case the eigenstates are given by Eq. (^) and the occupation 
numbers at a given energy are found from Eq. (|l0|). The key quantity in calculating n s (E) is the mean-squared 



eigenstate component 



(i) 



When V is strong enough and the energy £?W is not too close to the ground state, 



C 



(i) 



represents the spreading of the eigenstate over a large number of basis states. It is proportional to the strength 



function (introduced by Wigner |13| and also known as the local density of states), 



C, 



(i) 



p(E) 



(16) 



where p(E) is the eigenvalue density. The last equality in ( |16|) implies that some averaging over the the energy interval 
greater than the level spacing has been performed at E^ w E. It follows from numerical calculations Jl],||,fl| as well 
as from analytical considerations jlO 13 lj] that p w (E,k) is a bell-shaped function centered at E sw Ek- Near its 



maximum it depends only on the difference E — Ek, and can be described by the Breit- Wigner formula 



Pw{E, k) 



T/27T 



(E-E k ) 2 + T 2 /4 



(17) 
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The spreading width T is usually given by T ~ 2ir\(k'\V\k)\ 2 p(E). Therefore, the basis states are strongly mixed 
together by the residual interaction only locally, within the energy range T. 

The notions of the strength function and the spreading width become meaningful if the interaction is strong 
enough, and T ^ D, where D = l/p(E) is the mean level spacing (accordingly, |(fc'|U|/c)| 2 ^ D 2 ). This means that 
the number of basis states participating in a given eigenstate is large, N c ~ T/D ^> 1, or vice versa, a given basis state 
k contributes to a large number of nearby eigenstates with energies \E — E k \ ~ T. Apart from the smooth variation 

the statistics of the eigenstate components Ci is close to Gaussian ffl. In this situation it is appropriate 



of 



to call the stationary states of the system chaotic or compound eigenstates. 

Using Eq. (Jig) we can re-write expression din) for the occupation numbers in the integral form 



n.(E)= C£> (k\n s \k)p(E k )dE k n p w {E,k)(k\n s \k)dE k , (18) 



where we used the fact that p(Ek) ~ p(E) near the maximum of p w (E, k). The above representation is very similar to 
Eq. ( |l2|) for n s (T). Equation (|l8| ) describes averaging over the compound state strength function p w (E, k) of width T 
centered at energy E, whereas Eq. ( |l2| ) refers to a thermal average with the WT(E k ) probability density, whichpeaks 
near E(T). Of course, the width of the distribution WT{E k ) depends on the temperature [~ aT 3 / 2 for Eqs. (|) and 
(|l5|)], but if we choose the temperature by setting E — E(T), the two averages n s (E) and n s (T) should be close to 
each other, provided the widths of p w {E, k) and wriEk) are much greater than the multiparticle level spacing, and 
the difference between these widths does not exceed the single-particle energy interval in the system Q . 

In the next section we calculate thermally averaged occupation numbers and establish a relation between the 
effective temperature T and the excitation energy for the Ce atom. Numerical calculations will confirm that a 
temperature-based statistical theory agrees with the dynamic calculation, and describes well the peculiar behaviour 
of the occupation numbers in Ce. 

V. CALCULATION OF THERMALLY AVERAGED OCCUPATION NUMBERS 

A. Statistical model 

Let us now perform a statistical calculation of the occupation numbers for a system of n particles distributed over 
r orbitals with energies e s and degeneracies g s (s = 1, . . . , r). We will assume that the two-body interaction of any 
two particles in the orbitals s and p is U sp , where both the direct and exchange terms are included: 

Usp = al — - X ) ^ ^2^ Va > 3 f 3a ~ Va P a P) ■ ( 19 ) 

9s\9p sp) QeSj g gp 

Thus, U sp is averaged over the degenerate single-particle states within the orbitals s and p. 
The energy of a particular many-particle state k is now given by 



E k =^ s + ^^ j f r SP ' U SP , (20) 

s=l s=l p=s F 



where J\f s is an integer number of particles in the orbital s (0 < J\f s < g 8 ), and ^ s J\f s — n. The state k is specified by 

the orbital occupation numbers N s and is Gfe-degenerate, where G k = 111=1 (m, )• ^ course, Eq. (20) corresponds 

to the "diagonal" approximation [E k w H k k), since we neglect the effect of mixing of states due to the residual 
interaction (off-diagonal part of the Hamiltonian) . In the low-energy part of the spectrum this interaction pushes the 
eigenstates down with respect to their diagonal- approximation values because of the usual level-repulsion effect. 

In this model we cannot keep hold of the total angular momentum, so our calculation will yield quantities averaged 
over various angular momenta J and projections J z . However, it is relatively easy to ensure the conservation of parity. 
If the orbital s has a parity of P s (either 1 or —1), the parity of the multiparticle state k is Y\ r s =i ■ Therefore, one 
can easily select multiparticle states with a given parity when calculating statistical sums like that of Eq. (|lTJ) . 

In the canonical ensemble the probability of finding the system in the state k is given by 

w k = Z- 1 G k exp(-E k /T) , (21) 
where Z = ^ G k exp(-£; fe /T) , (22) 



G 



and the sum over k runs over all multiparticle states, possibly with a restriction on parity. The average occupation 
numbers n s (T) = AT s /g s are calculated as 



n s (T) = g; 1 J2^ k) w k , (23) 

fe 

where A/s is the number of particles in the orbital s in the multiparticle state k. In the diagonal approximation the 
energy of the system is related to the temperature as E = E(T), where 

£(T) =^£ feWfe , (24) 
fc 

is the canonical average. However, one can include the off-diagonal part of the Hamiltonian in the definition of 
temperature by introducing the energy shift Ae(T), 

E = E(T)-A E (T). (25) 

A e (T) is positive in the lower half of the spectrum, which means that the statistical effects of interaction between 
particles increase the effective temperature. At high temperatures Ae(T) = 2 < Hu — £w >^ where Hu — is the 
simple energy shift due to the non-diagonal matrix elements of the Hamiltonian H^ . This estimate is based on the 
mean energy of the components \k) in the eigenstate (Ekh = J2k Hkk\Ck | 2 = E^ + Ae (see Q). 

In general, the occupatio n nu mbers obtained from Eq. (E3) will be different from those predicted by the Fermi- 
Dirac distribution (see Sec. VB ). In appendix |a] we look at how the FDD (0) is derived from the canonical statistical 



distribution, Eq. ( J23| ) , and see what are the limitations on the interaction between the particles for the derivation to 
be valid. 



B. Numerical calculations 



To perform numerical calculations of the occupation numbers for Ce in the statistical model outlined above we use 



the same set of 14 relativistic orbitals as in the CI calculation described in Sec. Ill . The orbital energies are obtained 
as e s = (s\H c \s), where H c is the frozen Dirac-Fock Hamiltonian of the Ce 4+ core, and the averaged Coulomb matrix 
elements U sp are found from Eq. (|l9|). Their numerical values for the 7 lowest orbitals 4/ 5 / 2 , 4/ 7 / 2 , 6si/ 2 , 5d 3 / 2 , 
5d 5 / 2 , Qpi/2, and 6p 3 / 2 are given in Table Q. For excitation energies below the ionization threshold these orbitals are 
the most important. 

Using the statistical model formulae we have calculated the occupation numbers n s (T) (|23| ) and the energy of the 
system ( p4| ) as functions of T (see Figs. 3a and 3b). The relation between the energy of the system and temperature 
E = E(T) — Ae(T) can be inverted and used to plot the dependence of the statistical model occupation numbers 
as functions of the energy E. In this work we are mostly interested in the low-energy part of the spectrum, and we 
put Ae{T) — const to fit the true ground-state energy of the system at T = 0. In Fig. ^ we compare the results 
of the statistical model with the energy-averaged occupation numbers obtained from the CI calculation of the Ce 
eigenvalues and eigenstates, Eqs. (|J) and ([To|). We see that the complicated non-Fermi-Dirac energy dependence of 
the occupation numbers in Ce is reproduced well by the statistical model. 

To study the effect of the off-diagonal matrix elements of the Hamiltonian on the temperature we calculate the 
canonically-averaged mean energy of the system (|24|) using the set of the basis-state energies E k = -Hfcfc and that of 
the exact eigenstates E^\ The difference between these two mean energies is plotted in Fig. M. It is almost constant 
at small temperatures and follows 



A E (T) ~ ^ k \ (26) 

at large T (see §). 

Note that the energy shift Ae(T) is larger than the simple difference between the diagonal matrix elements and 
exact eigenvalues Hu — E^ . This is because the true occupation numbers ( [l0| ) even in the exact ground state are not 
integer (see Fig. |j) due to the admixture of higher configurations. This means that the effective temperature of the 
ground state is already not zero (see discussion in H). 
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VI. DISCUSSION AND CONCLUSIONS. 



Most importantly, the agreement observed in Fig. |] means that the interaction between the particles indeed 
introduces some kind of equilibrium in the system ("micro-canonical" distribution). Moreover, averaging over it 
yields results close to those over a canonical ensemble ([H]) , with the temperature chosen to reproduce the total energy 
of the system. This equivalence is always true for large system where any albeit weak interaction between particles 
leads to equilibrium. However, in a few-particle system the residual two-body interaction must be strong enough to 
produce chaotic eigenstates and facilitate statistical description (see [[15 16 where criteria for the interaction strength 
are discussed). 

We would like to reiterate, that although the temperature-based description is apparently applicable to our 4- 
particle system, the orbital occupancies could not be described by the FDD (Fig. ||). The FDD is inapplicable to 
our system because of the strong interaction between the particles [second term on the right-hand side of Eq. (^p|)1. 
However, the deviation from the FDD is determined not by the magnitude of U sp , but rather by the size of their 
fluctuations. To see this assume for a moment that all two-body matrix elements are the same, U sp = U . In this 
case the double sum in Eq. ( pp| ) just shifts all energies by ^N(N — 1), and the statistical properties of the system 
are the same as for noninteracting particles. If U sp are different for different orbital pairs sp one can still introduce 
some average interaction U and subtract this "background" interaction from the interaction term in Eq. (pfj| ). This 
procedure effectively suppresses the interaction term, since the summands in expressions like ^2 s<p (U sp — U)J\f p have 
different signs. Note that the introduction of (energy-dependent) U is equivalent to a local mean-field approximation 
(see also Appendix ^|). This approximation is good if the fluctuations of U sp from one orbital to another are relatively 
small. 

An instructive example was provided in [^|Q| . In these works a model of random two-body interactions was explored 
numerically and a good Fermi-Dirac-like behaviour of the occupation numbers was observed for as little as 4 particles 
distributed among 11 orbitals. However, this regime was achieved for the relatively small two-body matrix elements 
with mean zero and r.m.s.V ~ O.ldi, where d\ is the mean level spacing between the single-particle orbitals. On the 
other hand, for smaller two-body interaction strengths the occupation numbers distribution was not smooth, and did 
not agree well with the Fermi-Dirac formula, because the statistical equilibrium needed was not achieved. 

The situation in the Ce atom is different. The 4/ orbital has a much smaller radius than any other orbital, and the 
Coulomb interactions have a hierarchy of scales: 

Uifif > Uif s > U sp , (27) 

where s and p are orbitals other than 4/. Indeed, the Coulomb interaction between the electrons is determined 
mostly by the mean radius of the largest orbital. Thus, for the two orbitals s and p such that r s < r p , the Coulomb 
interaction is U sp « e 2 /r p (this formula is exact for the two electrons distributed uniformly over the surfaces of two 
spheres of radii r s and r p ). Because of ( p7j ) the interaction term in Eq. ( |20| ) fluctuates strongly with the change of 
the occupation number of the 4/ orbital, and in effect there is no good mean-field approximation for the excitation 
spect rum of Ce. The quasiparticle orbital energies e s can still be obtained in the statistical model for Ce by means of 
Eq. ( AlC ). They are plotted in Fig. ^c and show considerable variation with temperature. However, even when we 



use these energies instead of the Hartree-Fock ones for plotting the occupation numbers, the resemblance to the true 
FDD is only marginally better than that in Fig. ||. 

The absence of reasonably defined quasiparticle orbitals and the ensuing distortion of the Fermi-Dirac distribution 
are features outside the usual Migdal theory of normal finite Fermi systems (TFFS) Jl7j. Its breakdown in the Ce 
atom can be associated with the open-shell structure of the atom (nearly degenerate "ground state" ) and a clear lack 
of symmetry of the ground state with one removed particle. As a result, single-particle excitations above the ground 
state do not carry good quantum numbers (like the momentum in an infinite Fermi system, or the angular momentum 
in a spherically symmetric finite system). Moreover, even at low energies (few eV) the single-particle excitations have 
large widths associated with their decay into multiply excited configurations (the spreading width V ~ 2 eV [Q), 
largely because such decay is not really limited by any selection rules (only the trivial total angular momentum and 
parity are conserved). 

It has been proposed in [|| that for finite Fermi systems similar to the Ce atom, characterized by the dense spectra 
of chaotic multiparticle eigenstates, a statistical theory alternative to the standard TFFS can be developed based on 
the properties of these eigenstates. Most importantly, the existence of the chaotic eigenstates and the equilibrium 
this introduces in the system is ensured by the sufficiently strong interaction between particles. 

This concept of the interaction-driven equilibrium is supported by our present results. We have shown that this 
equilibrium can be described in terms of usual statistical parameters, such as the temperature, even though some of 
the system's properties are very different from those usually expected in Fermi systems. For example, the statistics 
of the occupation numbers cannot be described by the Fermi-Dirac formula. 
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APPENDIX A: 



DISTRIBUTION OF OCCUPATION NUMBERS FOR A CANONICAL ENSEMBLE OF 
FINITE SYSTEMS OF INTERACTING FERMIONS 



Consider a quantum system which consists of a number of single-particle discrete states a (we will also call them 
orbitals here) with energies e a (a = 1, . . . ,m) filled with n < m Fermi particles. The multiparticle states k of the 
system arc identified by specifying the occupation numbers n a = or 1 of the orbital a, n a = n - The total 
number of multiparticle states N — (™) is quite large, even for moderate m and n. If we allow for the interaction 
between the particles, the energy of the state k is given by 



E k = s n n 



ol(3 



(Al) 



where U a fj includes both the direct and exchange interaction between the particles in a and j3, and U aa = 0. 

Of course, the states k are not eigenstates of the system. However, if the interaction between the particles is not 
too strong, we can use them for averaging over the canonical ensemble with probabilities w k = Z^ 1 exp(—Ek/T), 
where Z = J2 k ex P( — E k /T) is the partition function, and the sum runs over all N multiparticle states pq . It will 
be convenient to show explicitly that the partition function depends on the energies of the orbitals e%, . . . ,e m = {e}, 
number of particles n, and temperature T: Z = Z({e},n,T) (and, strictly speaking, on the interactions U a p, if they 
are not zero). 

Using the canonical probabilities we can express the mean occupancy of the orbital a as 



exp(- J B fe /T) 



Z k eM-E k /T) 
E fc( «) *M-Ek/T) 



£*(„) eM-Ek/T) + J2 k(a) eM-Ek/T) ' 



(A2) 
(A3) 



fe(a) 



where = 1 or 0, depending on whether a is occupied or empty in the state k, and we split the sum over k into 
two sums over the states k(a) where a is occupied, and k(a), where it is empty. 

Let us first consider the case of noninteracting particles (U a p = 0). It is easy to see that the first sum is then equal 
to Z{{e}' a , n — 1, T) exp(— e a /T) and the second one is Z({e}' a ,n, T), where {e}^ is the set of m — 1 orbitals obtained 
by discarding a from {e}. Equation can then be written as 



Z({e}' a ,n,T) 



Z{{e} c 



1,T) 



exp(e Q /T) 



This equation is very similar to the Fermi-Dirac formula (|l|), if we introduce the chemical potential /i by 



Z({e}' a ,n,T) 
Z({e}' a ,n-1,T) ~ 



exp(-/x/T) 



(A4) 



(A5) 



The problem is that this ratio on the left-hand side in fact depends on which orbital a is deleted from the set {e} to 
form {e}' a , and so does the "chemical potential" fj,. If we write Eq. (AA ) for np with (3 ^ a, the set {e}^ will produce 
a different ratio Z({e}'p, n, T)/Z({s}'p,n— 1, T), and as a result, a different value of fj,. However, {e}^ can be obtained 
from {e}„ by simply moving the orbital energy from ep to s a . S o, the difference between the values of \x for different 



orbitals can be probed by calculating the derivative of Eq. (A5) with resp ect to the ene rgy of some orbital (3. Using 
the relation Tip — —TdZ/dep, valid for noninteracting particles (see Eqs. (Al) and (|A2{)) we obtain 



9/i 



Jn-l) 



(A6) 
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where ' and 



(n-l) 



are the mean occupation numbers for n and n — 1 particles distributed among the {e}^ orbital 



set. It is obvious that the right-hand side of Eq. (A6) is larger near the Fermi level, \eg — /i| < T, and is almost 



zero outside this interval. We can estimate that the total difference between the value of n for orbitals well below the 
Fermi level and those well above it is 



5n = J (7 



in) 



.(n-l) 



den — di 



,(n) 
l /3 



_(n-l) 



^wd 1 [n-(n-l)]=d 1 
«i 



(A7) 



where di is the mean spacing between the single-particle orbitals. Thus, in a system with discrete orbital energies the 
chemical potential can be considered as constant to within ~ d\ accuracy. At finite temperatures the width of the 
smoothened Fermi-Dirac "step" is of the order of T, therefore \i = const is valid for T 3> d\ (or for /1 ^> d\). Note 
that this condition means that the number of orbitals within the Fermi-Dirac "step" is large. 

For the interacting particles the sum J2k(a) ex P(~-E'fe/7 1 ) which gives rise to exp(e Q /T) in Eq. (A4) can be written 



eM-Ek/T) = ex P 



k(a) 



k{a) 



x exp 




Uapnp 



(A8) 



In this formula the last exponent contains the energy of the particle in a, which depends on the occupancies ng of 
the other orbitals. When summation over k(a) is carried out the exponent is averaged over different distributions of 
n — 1 particles among all orbitals but a. The result can be presented approximately as 



Z({e} a ,n-l,T)exp(-e Q /T) , 
where we replaced the averaged exponent by the exponent containing the quasiparticle energy, 



(A9) 



(A10) 



and the mean values ng are, strictly speaking, different from th ose from Eq. ([A4), as one particle has always been 
kept in a in the sum (A8). Therefore, the "Fermi-Dirac" anzats (A4) holds for the interacting parti cles, if we replace 
the single-particle energies e a with the tem perature-de pen de nt qu asiparticle energies e a from Eq. ( AlC ) 

Note that the transformation of Eq. ( A8) into Eqs. ( A9), ( A10) is exact up to first order in Hp, and to all orders, if 
U a p = U for all orbitals. In the latter + (n — 1)U is a trivial redefinition of the single-particle energies. 

That is why replacing e a with the quasiparticle energies e a is a valid operation, unless the interactions U a g fluctuate 
strongly, and the number of active particles is small. In the latter case the mean-field approximation is inadequate 
and the introduction of quasiparticles is not very meaningful. 
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TABLE I. Single-particle energies and Coulomb matrix elements for the valence and lowest excited orbitals in Ce. 



Orbital 








Coulomb matrix elements U sp , a.u. 






nlj 


a.u. 


4/5/2 


4/7/2 


6S1/2 


4^3/2 


4d 5/2 


6^1/2 


6P3/2 


4/s/2 


-1.564 


0.791 


0.800 


0.260 


0.423 


0.422 


0.223 


0.216 


4/7/2 


-1.551 


0.800 


0.787 


0.259 


0.428 


0.416 


0.223 


0.215 


6S1/2 


-0.876 


0.260 


0.259 


0.199 


0.231 


0.230 


0.162 


0.158 


4^3/2 


-1.156 


0.423 


0.428 


0.231 


0.330 


0.331 


0.200 


0.198 


4d 5 /2 


-1.141 


0.422 


0.416 


0.230 


0.331 


0.325 


0.204 


0.195 


6P1/2 


-0.714 


0.223 


0.223 


0.162 


0.200 


0.204 


0.169 


0.157 


6P3/2 


-0.691 


0.216 


0.215 


0.158 


0.198 


0.195 


0.157 


0.156 



FIG. 1. Eigenvalue densities for the even states of Ce averaged over the energy interval AE = 0.05 a.u. a. The upper solid 
curve is for J = 4, and lower curves correspond to successively increasing values of J. Note that all densities have similar shapes. 
Dotted curve is the analytical fit for J — 4: p(E) = poexp \a{E — -Eg) 1 ^ 2 ] , where po = 119, a = 12.3 a.u. and E g — —2.91 a.u. 
b. Total level density of the even states and fit with po = 27, a = 13.0 a.u., and E g — —2.91 a.u. 



FIG. 2. Occupation numbers n s (E) (see Eq. (|To|)), for the even states of Ce at the excitation energy of E — E g — 4.5 eV 
versus the single-particle energies e s of the orbitals (a), and quasiparticle energies e s (b). Diamonds connected by the dashed 
line represent the result of our statistical calculation with T = 0.097 a.u. 



FIG. 3. Temperature dependence of the occupation numbers for the orbitals g a n s (T) f a), energy of the system E(T) (b) 
and quasiparticle energies i s (c), calculated in the statistical model of the Ce atom, Eqs. (bo|)-(E4). The occupation numbers 
shown for the 4/, 5d and 6p subshells are sums of those of their fine-structure sublevels j — I ± 



FIG. 4. Comparison of the orbital occupancies g 3 n a (E) obtained from the exact diagonalization (solid and dash-dot lines) 
with g a n s (T(E)) obtained from our statistical theory (dotted lines). 



FIG. 5. Difference between the mean values of the energy obtained from the canonical distribution using the energies of the 
basis-states, Et = Hkk, and those of the exact eigenstates. The dash-dot line is the high-temperature analytical approximation 
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Figure 2. 
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Figure 5. 



